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Numerical  Simulation  of  Boundary  Layer  Transition* 

by 

Steven  A.  Orszag 

Consultant,  Flow  Research  Company,  Kent,  Washington  98031 


Abstract 

The  purpose  of  this  research  was  to  develop  numerical  simulations 
of  boundary  layer  transition.  One  of  the  essential  aspects  of  the  work 
was  to  investigate  carefully  the  appropriate  boundary  conditions  and 
mathematical  models  that  should  be  employed  in  the  study  of  transition 
processes.  The  effects  of  three  dimensionality,  spanwise  wavelength 
selection,  and  streamline  curvature  were  also  considered. 
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1 . Introduction 

In  this  report,  we  summarize  the  results  obtained  under  Office  of 
Naval  Research  Contract  No.  N00014-76-C-0263 , Task  No.  NR  061-233,  ARPA 
Order  No.  2924.  The  purpose  of  this  research  was  to  develop  numerical 
simulations  of  boundary  layer  transition.  One  of  the  essential  aspects 
of  the  work  was  to  investigate  carefully  the  appropriate  boundary 
conditions  and  mathematical  models  that  should  be  employed  in  the  study 
of  transition  processes.  The  effects  of  three  dimensionality,  spanwise 
wavelength  selection,  and  streamline  curvature  were  also  considered. 

The  direct  solution  of  the  Navier-Stokes  equations  for  the  study  of 
transition  involves  much  computation  because  three-dimensional  effects 
are  important  (see  section  3),  and  high  resolution  is  required  in  the 
downstream  and  boundary  layer  directions.  To  economize  in  the  solution 
of  transition  problems,  we  have  developed  a nonlinear  stability  theory 
that  seems  to  account  for  the  principal  effects  observed  in  our  direct 
numerical  simulations.  This  Galerkin  approximation  is  discussed  in 
section  4. 

In  this  report  we  briefly  summarize  our  results.  Detailec  exposi- 
tions of  our  results  are  given  in  the  publications  cited  in  the  references. 
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2.  Boundary  Conditions 

The  basic  flow  geometry  of  the  numerical  simulations  of  transition 
reported  in  section  3 is  shown  in  figure  1.  The  flow  enters  the  com- 
putational domain  at  x = 0 and  exits  at  x = L (normally  L = 50  cm  ) . 

The  free  stream  is  at  z = 00  , while  the  rigid  wall  is  at  z = 0 . In 
the  spanwise-y  direction,  the  unperturbed  boundary  layer  flow  is  assumed 
uniform  while  the  perturbed  (transitional)  flow  is  assumed  to  satisfy 
periodic  boundary  conditions.  A typical  calculation  reported  in  section 
3 involves  the  use  of  257  grid  planes  to  resolve  the  downstream  x direction, 
8 Fourier  modes  to  resolve  the  spanwise-y  direction,  and  33  Chebyshev 
polynomials  to  resolve  the  boundary-layer  z direction.  We  discuss  the 
boundary  conditions  applied  in  each  of  the  coordinate  directions  in  more 
detail  below. 

Spanwise  Direction  - Periodicity 

The  assumption  of  periodicity  in  the  spanwise  direction  is  con- 
sistent with  the  Navier-Stokes  equations  in  the  sense  that  if  the 
perturbed  boundary  layer  flow  is  initially  periodic  in  v at  t = 0 , 
then  the  flow  stays  periodic  in  y with  the  same  period  for  all  later 
times.  This  periodicity  does  not  prevent  the  formation  of  localized 
vortex  structures  in  the  flow  during  the  process  of  transition — localized 
vortex  structures  that  do  appear  during  time  evolution  must,  however, 
appear  periodically  in  the  y direction.  The  assumption  of  spanwise 
periodicity  does  not  seem  to  cause  any  serious  difficulty  or  physical 
inconsistency  in  comparing  our  computations  with  laboratory  experiments 
like  those  of  Klebanoff,  Tidstrom,  and  Sargent  (1962). 

Boundary  Layer  Direction  - Mapping 

We  have  discovered  that  an  algebraic  mapping  of  the  semi-infinite 
domain  above  the  flat  plate  located  at  z = 0 into  the  finite  interval 
(-1,1)  gives  remarkably  efficient  and  accurate  numerical  computations. 

The  algebraic  mapping  is 


z 


z - H 
z + H ’ 


(2.1) 


where  1!  is  a scale  paramete 
which  the  unperturbed  velocity 


aally  chosen 
is  one-half  the 


to  be  twice  the  height  at 
free-stream  value  (for 


« 


i 


sec 
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the  Blasius  boundary  layer  this  height  is  approximated  well  by  twice 
the  displacement  thickness  of  the  boundary  layer). 

The  mapping  (2.1)  transforms  the  interval  0 z < “ into  the 
finite  interval  -1  £ z < 1 . A nice  feature  of  the  mapping  (2.1)  is 
that  the  rule  for  transformation  of  derivatives  is  simply 

9u  _ (1  - z)~  3u  n 0. 

3z  2H  32  • u*“"' 

Basically,  the  mapping  (2.1)  is  just  an  analytical  device  to  introduce 
nonuniform  resolution  in  the  z direction.  If  a numerical  scheme  with 
uniform  resolution  in  z is  used,  too  much  information  is  wasted  in 
the  free  stream,  where  the  solution  changes  slowly.  This  puts  undue 
burdens  on  the  numerical  resolution  closer  to  the  wall,  where  the  action 
lies. 

Grosch  and  Orszag  (1977)  have  given  a detailed  account  of  our 
analysis  of  mapping  in  this  coordinate  direction.  They  give  a large 
variety  of  examples  and  explain  in  detail  why  the  z direction  is  ripe 
for  mapping  techniques.  Here  we  give  just  two  examples  that  illustrate 
how  well  mapping  works  in  the  z direction.  The  two  examples  are  the 
Orr-Sommerf eld  equation  for  the  linear  eigenmodes  of  Blasius  boundary 
layer  flow  and  the  Falkner-Skan  equation  for  the  stream  function  of 
boundary  layer  flow  on  wedges  (with  a pressure  gradient). 

In  table  I,  we  list  the  most  unstable  eigenvalue  of  Blasius  flow 
at  a Reynolds  number  R = Ux/v  = 580  , where  x is  the  distance 
from  the  leading  edge  of  the  flat  plate,  with  a (dimensionless)  wave  number 
of  0.179.  Two  numerical  schemes  are  compared:  a Chebyshev  polynomial 
solution  (Orszag,  1971)  of  the  Orr-Sommerf eld  equation  truncated  to  the 
domain  0 z £ H for  various  values  of  H with  the  boundary  condition 
f'  + 0.179f  = 0 applied  to  each  of  the  eigenfunctions  at  z = H ; and  a 
Chebyshev  polynomial  solution  of  the  same  problem  with  the  mapping  (2.1) 
with  scale  factor  H = 1 . Note  that,  according  to  the  discussion  following 
(2.1),  the  optimal  choice  for  H is  about  3.4,  so  we  could  anticipate 
even  greater  accuracy  by  such  a choice.  (Here  and  in  the  rest  of  our 
work  we  non-dimensionalize  z by  the  length  scale  vx/U  , where  U 
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is  the  free-stream  velocity,  so  the  displacement  thickness  of  the 
Blasius  layer  is  1.7.)  Also  note  that  it  is  not  necessary  to  apply  any 
boundary  conditions  at  all  at  z = 1 when  using  the  mapping  since  eq. 
(2.2)  ensures  that  all  nonsingular  derivatives  vanish  at  z = 1 after  , 
mapping . 

Table  I Eigenvalues  of  the  Orr-Sommerfeld  Equation  For  Blasius  Flow 


Mapping 


Restricted  Domain 
0 < z < H 


H N 

Number  of 
Scale  Chebyshev 

Factor  Polynomials 


C 

Eigenvalue 


20  44  0.360213  + iO. 006671 

30  44  0.364041  + iO. 008113 


Algebraic  Mapping  1 

z = (z-H) / (z+H)  1 

Exact  Eigenvalue  c 


26  0.364147  + iO. 008007 

34  0.364121  + iO. 007958 

0.364123  + iO. 007960 


In  Table  II  we  list  some  errors  in  the  value  of  u(l)  , which  is 
the  non-Jimensionalized  x component  of  the  velocity  at  a non-dimensionalized 
height  of  1 above  the  rigid  wall,  as  determined  by  numerical  solution  of 
the  Falkner-Skan  equation  with  only  11  grid  points.  Three  different 
pressure  gradients  are  listed  in  the  table.  For  each  case,  the  alge- 
braic mapping  achieves  much  better  accuracy  than  that  achieved  by 
simply  restricting  the  integration  domain  to  0 < z £ H . In  table  II, 
the  points  labelled  * indicate  that  no  acceptable  (monotonic)  solution 
to  the  discretized  Falkner-Skan  equation  exists  (with  only  11  grid  points). 

Unfortunately,  the  technique  of  mapping  works  well  only  for  problems 
where  the  solution  is  "simple"  at  infinity  (see  Grosch  and  Orszag  (1977) 
for  a variety  of  examples).  Thus,  boundary  layer  flows  are  "simple”  at 
infinity  because  perturbations  that  are  large  near  the  wall  die  away 
exponentially  fast  as  z increases.  For  example,  if  the  amplitude  of 
a linear  mode  of  the  Orr-Sommerfeld  equation  oscillates  as  elj,x  , then 

“OlZ 

its  amplitude  decreases  like  e for  large  z . On  the  other  hand, 
the  downstream  x direction  is  not  one  in  which  the  flow  becomes  "simple"; 
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rather,  we  expect  that  sufficiently  far  downstream  the  flow  will  become 
turbulent.  Mapping  the  semi-infinite  x domain  into  a finite  interval 
without  considering  the  outflow  boundary  conditions  will  not  work. 

Table  II  Errors  in  Solution  of  the  Falkner-Skan  Equation 


Mapping 

H 

3 

Error  in  f ' (: 

Restricted  Domain 

1.0 

0 

0.54 

2.5 

0 

0.052 

5.0 

0 

0.00013 

7.5 

0 

* 

Algebraic  Mapping 

1.0 

0 

0.0000058 

2.5 

0 

-0.00000047 

5.0 

0 

-0.0000028 

7.5 

0 

-0.0000040 

Restricted  Domain 

5.0 

-0.1 

0.00042 

7.5 

-0.1 

jV 

Algebraic  Mapping 

5.0 

-0.1 

-0.0000046 

7.5 

-0.1 

-0.0000064 

Restricted  Domain 

2.5 

0.1 

0.032 

5.0 

0. 1 

0.00014 

7.5 

0.1 

A 

Algebraic  Mapping 

2.5 

0. 1 

-0.0000055 

5.0 

0.1 

-0.000011 

7.5 

0.1 

-0.000016 

^ Solutions  of  the 

Falkner-Skan  equation 

obtained 

by  using  11  equally 

spaced  grid  points 

in  the  appropriate  coordinate 

system.  Here  * 

indicates  that  no 

acceptable  (monotonic) 

solution 

was  found.  Also, 

f'(l)  = u(l)/U  , the  relative  x velocity  at  a non-dimensional  height 
of  1;  770  is  the  included  angle  of  the  analogous  wedge  flow. 


Inflow-Outflow  Boundary  Conditions 

Orszag  (1974)  has  shown  that  the  viscous  Navier-Stokes  equations  are 
well  posed  if  all  three  velocity  components  are  specified  at  inflow  and 
outflow  boundaries.  Orszag  also  showed  that  the  inviscid  Navier-Stokes 
equations  are  well  posed  if  all  three  velocity  components  are  specified 
at  inflow  points  and  if  either  the  normal  component  of  the  velocity  or 
the  pressure  is  specified  at  an  outflow  point.  We  have  investigated  in 


Flow  Research  Report  No.  80 
May  1977 


detail  the  effects  of  various  boundary  conditions  on  solutions  of  the 
Navier-Stokes  equations  and  have  developed  new  boundary  conditions  that 
minimize  the  effect  of  arbitrarily  specified  outflow  boundary  conditions 
on  the  flow  in  the  interior  of  the  domain  of  interest.  A detailed 
report  on  our  studies  is  being  prepared  (Orszag  and  Israeli,  1977);  a 
preliminary  outline  of  our  results  was  published  by  Orszag  (1976). 

A simple  model  of  viscous  effects  is  given  by  the  linearized 
Burgers'  equation: 

lH+u3H=vl!ii  (2  ' 

3t  3x  V „ 2 ’ 

ax 

where  U is  a constant  advecting  velocity  and  V is  the  kinematic 
viscosity.  Suppose  we  wish  to  solve  eq.  (2.3)  on  the  semi-finite 
domain  0 _<  x < °°  , and  suppose  that  U > 0 so  that  x =0  is  an 
inflow  point.  Let  us  see  what  happens  if  the  domain  is  truncated  to 
the  finite  interval  0 £ x 1 . If  v is  very  small,  the  Laplace 
transform  solution  of  the  initial-value  problem  for  eq.  (2.3)  is  composed 
of  modes  of  the  form 

\v4-fTt 

u(x,t)  = e X , (2.4 


where  the  dispersion  relation 


vX  - XU  - a = 0 


determines  as  a function  of  c . Thus,  if  V«  1 , there  are  two 

possible  modes  given  approximately  by 


The  first  mode  X^  is  a viscous  mode  that  decays  as  -x  increases.  The 
second  mode  X£  is  a mode  that  describes  propagation  from  x = 0 to 


x “ + » . 
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The  general  solution  to  eq.  (2.3)  is  of  the  form 

A _ x X„x 

at  1 ° 

u(x,t)  = dae  A^(o)e  + A?(o)e  . (2.8) 

When  the  semi-infinite  domain  0 x < °°  is  truncated  to  the 
finite  interval  0 £ x <_  1 , it  is  necessary  to  impose  boundary  conditions 

at  both  x = 0 and  x = 1 , say  u(0,t)  = f(t)  and  u(l,t)  = 0 . Since 

X.  is  large  and  positive,  it  follows  from  eq . (2.8)  that  the  effect 

of  this  mode  at  x = 0 is  exponentially  smaller  than  its  effect  at 

x = 1 . Therefore,  neglecting  the  mode  X , we  find  from  eq.  (2.8) 
that  the  boundary  condition  u(0,t)  = f(t)  determines  A„,(a)  . On 
the  other  hand,  the  boundary  condition  u(l,t)  = 0 determines  A^(c) 
in  terms  of  A0(a)  : 

A1  = _A2  exp( ' ■ “ V • (2,9) 

Now  we  can  compare  the  finite  an!  semi-infinite  problems.  In  the 
semi-infinite  problem  the  mode  X.  does  not  appear  at  all,  but  the 
solution  is  otherwise  the  same  as  for  the  finite  problem.  Therefore, 
according  to  eqs.  (2.8)  and  (2.9)  the  effect  of  the  boundarv  condition 
imposed  at  x = 1 is  confined  to  a boundarv  layer  of  width  of  order 
v/U  near  x = 1 . 

The  boundary  condition  u(l,t)  = 0 affects  the  solution  only  in 
the  thin  boundary  layer  with  width  of  order  v/U  . It  is  possible  to 
decrease  the  effect  of  the  artificial  boundary  at  x = 1 still  further 
by  imposing  a boundary  condition  of  the  form 

£— r u ( 1 , t ) = 0 (2.10) 

a K 

ox 


at  x = 1 . In  this  case,  the  solution  for  A^  that  corresponds  to 
eq.  (2.9)  is 
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straightforward  application  of  u(l,t)  = 0 . The  only  danger  in  applying 
higher  order  conditions  like  eq.  (2.10)  with  k large  is  that  numerical 
instability  may  be  induced  by  high  order  extrapolations  near  the  outflow 
boundary.  With  numerically  stable  schemes,  however,  higher  order  boundary 
conditions  should  be  more  accurate  than  lower  order  boundary  conditions. 

We  can  make  a similar  analysis  of  a model  problem  that  bears  much 
closer  relation  to  incompressible  fluid  mechanics  than  does  eq.  (2.3). 

The  model  includes  the  effect  of  pressure,  which  is  essential  for  under- 
standing the  nature  of  the  Navier-Stokes  equations.  The  model  consists 
of  the  two-dimensional  linearized  Navier-Stokes  equations  with  a solution 
assumed  to  be  of  the  special  form 

u = u(x,t)elky  , v = v(x,t)elky  . (2.12) 

The  model  equations  are  given  by 

u + Uu  + ikVu  = - p + v(u  - k^u)  , (2.13) 

t x rx  xx 

v + Uv  + ikVv  = - ikp  + v(v  - k'v)  , (2.14) 

t x r XX 

ux  + ikv  = 0 , (2.15) 

where  we  assume  that  U > 0 so  x = 0 is  an  inflow  boundary  and  V 
is  arbitrary.  We  want  to  study  the  effect  of  truncating  the  semi- 
infinite domain  0 < x < 00  into  the  finite  domain  0 <_  x <_  1 . To  do 
this,  we  analyze  the  Laplace  transform  solution  of  the  initial  value 
problem  for  eqs.  (2. 13)-(2. 15) . 

The  Laplace  transform  solution  of  eqs.  (2. 13)— (2. 15)  is  composed  of 
modes  of  the  form 


, . Ax+ot 

u(x,t;  = ue  , 

(2.16) 

, Ax+ot 

v(x,t)  — ve  , 

(2.17) 

r 
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where  the  dispersion  relation 


(X2  - k2)(o  + XU  + ikV  - vX2  + vk2)  = 0 


(2.18) 


determines  X as  a function  of  O . The  four  roots  of  eq.  (2.18)  for 
given  a , k , U , and  V are 


X.^  ~Fk.,  X ^ — k , 


and  \_  , X,  are  the  solutions  of  the  quadratic  equation 

J *4 


vX  - XU  - (a  - ikV  - vk“)  = 0 

When  v«l  , the  asymptotic  behaviors  of  X^  and  A ^ are  given  by 

o+ikV 


(2.19) 


X,  - U/v  , X, 

J *4 


u 


(2.20) 


The  general  solution  of  eqs.  (2. l3)-(2. 15)  is  of  the  form 

'•  x X,x  \ x X x 

u(x,t)  = dc  e'  A,  (o)e  ^ + A.,(a)e  “ + A^(')e  + A.  (c)e  "*  . (2.21) 


In  the  semi- inf inite  domain  only  modes  2 and  4 survive  if  the 
solution  is  bounded  at  « , so  eq.  (2.21)  reduces  to 

X.x  X,x 

n t 2 -* 

u(x,t)  = da  e A-(a)e  + A • (a)e  . (2.22) 

z. 

In  the  finite  domain  0 < x £ 1 , all  four  modes  are  present  to  some 
extent  and  are  determined  by  the  boundary  conditions.  The  contribution 
of  mode  3 decays  in  a viscous  boundary  layer  around  the  outflow  boundary 
at  x = 1 (as  for  the  linearized  Burgers'  equation),  so  if  v is 
small,  this  mode  causes  negligible  error  throughout  the  interior  of  the 
interval  0 <_  x < 1 . On  the  other  hand,  the  mode  X ^ = +k  is  a pressure 
mode  that  decays  appreciably  only  over  the  distance  1/k  . so  if  k is 
not  large,  its  effect  may  be  large  over  the  whole  computational  domain. 

It  is  mode  1 that  is  the  origin  of  the  principal  differences  between  the 
semi-infinite  flow  case  and  the  finite-domain  simulation. 
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There  are  several  ways  to  minimize  the  effect  of  mode  1.  One  way  is  to 
assume  that  the  computational  domain  has  a large  aspect  ratio  L/H  >>  1 , 
where  L is  the  length  of  the  computational  region  and  H is  its 
height.  In  this  case  the  smallest  allowable  transverse  wavenumber 
k satisfies  k + 2tv/H  » 1/L  , so  that  mode  1 does  in  fact  "boundary 
layer"  (in  a distance  of  order  H ) about  the  outflow  boundary  x = L . 

Thus,  so  long  as  L - x » H , the  effect  of  arbitrarily  specified 
outflow  boundary  conditions  at  x = L is  negligible. 

Another  way  to  limit  the  effect  of  the  outflow  boundary  conditions 
is  to  modify  the  equations  of  motion  in  the  neighborhood  of  the  outflow 
boundary  so  that  the  "bad"  mode  1 boundary  layers  in  a thin  boundary 
layer  near  x = 1 . We  therefore  make  the  wavenumber  of  mode  1 become 
A^  » 1 . Three  ways  to  modify  the  Navier-Stokes  equations  to  accomplish 
this  wavenumber  are: 

(i)  Replace  the  incompressibility  condition  u + v =0  by  u 

■+  x y x 

+ Vy  = y v*n  , where  y»l  , and  n is  the  outward  normal. 

(ii)  Modify  the  pressure  term  in  the  Navier-Stokes  equations  to  be 
-^p  + ynp  . 

(iii)  Introduce  pseudo-compressibility  terms  in  the  Navier-Stokes 
equations  near  the  outflow  boundary. 

For  example,  if  we  replace  eq.  (2.15)  by 

u + ikv  = yu  , (2.23) 

x 

modes  3 and  4 of  eqs.  (2.13),  (2.14),  and  (2.23)  are  unchanged  from  eq. 
(2.20),  while  modes  1 and  2 are  replaced  by 

~ V , A2  ~ - k2/u  , (2.24) 

when  y » 1 . Thus,  the  size  of  A^  is  increased,  and  mode  1 boundary 
layers  near  x = 1 . To  minimize  the  effect  of  the  outflow  boundary  at 
x = 1 , we  choose  y(x)  such  that  the  effects  of  mode  1 are  ocalized 
x = 1 . 


near 


▼ 
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This  technique  has  been  tested  on  the  problem  of  low  Reynolds 
number  flow  past  a cylinder  with  mixed  success.  At  a Reynolds  number 
of  four,  when  a truncation  of  the  computational  domain  to  1 r £ r^ 
is  used,  the  results  for  the  drag  coefficient  are  as  follows.  Using 

the  Imai  asymptotic  behavior  of  the  vorticity  and  stream  function  as 
boundary  conditions  at  r^  , we  obtain  = 4.435  if  r^  = 21.5  , 

CD  = 4.441  if  r^  = 35.0  , and  CD  = 4.442  if  r^  = 59.1  . Using 
free-stream  boundary  conditions  instead,  we  obtain  = 5.13  if 
= 21.5  , CD  = 4.66  if  r^  = 59.1  , showing  the  very  large  errors 
that  can  result  from  use  of  free-stream  boundary  conditions.  On  the 
other  hand,  the  use  of  a modified  Navier-Stokes  equation  of  the  form  (i) 
with  derivative  boundary  conditons  applied  at  r and  a quartic  power 
law  for  u(r)  beyond  0.8r  , we  obtain  CL  = 4.465  if  r = 21.5  . 

However,  if  we  use  non-dif f erentiated  boundary  conditions  on  the  outer 
radius  r , then  the  best  choice  of  p(r)  gives  Cn  =4.73  at  r = 21. 
Thus,  with  differentiated  boundary  conditions,  the  new  methods  appear  to 
give  superb  results,  while  with  non-dif ferentiated  boundary  conditions, 
the  results  are  less  spectacular.  A full  discussion  of  these  results 
and  methods  is  given  by  Orszag  and  Israeli  (1977). 
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3 . Numerical  Simulations  of  Boundary  Layer  Transition 

In  this  section  we  describe  some  numerical  results  for  boundary 
layer  transition.  In  our  earlier  work  on  transition  (Orszag,  1974), 
we  solved  a set  of  "parabolized"  Navier-Stokes  equations  in  which  the 
x derivatives  of  the  viscous  dissipation  terms  were  dropped  on  the  basis 
that  they  were  expected  to  be  roughly  100  times  smaller  than  the 
z-derivative  contribution  to  the  viscous  terms  in  linear  Tollmien- 
Schlicting  waves.  The  advantage  of  the  parabolized  Navier-Stokes 
equations  is  that  they  allow  the  use  of  inviscid  boundary  conditions  at 
the  outflow  boundary  and  prevent  the  occurrence  of  thin  boundary  layers 
near  the  outflow  boundary.  In  the  present  work,  we  have  compared  the 
results  of  the  parabolized  Navier-Stokes  equations  with  numerical  solutions 
of  the  full  Navier-Stokes  equations  and  have  found  good  agreement  out- 
side a thin  boundary  layer  near  the  outflow  boundary.  The  results 
presented  below  were  all  obtained  by  solution  of  the  full  Navier-Stokes 
equations . 

Numerical  Methods 

The  Navier-Stokes  equations  were  solved  by  using  257  staggered 
grid  planes  in  x with  a compact  mesh  fourth-order  difference  operator 
(Orszag  and  Israeli,  1974).  Inflow  boundary  conditions  at  x = 0 were 
that  all  three  velocity  components  be  specified  (see  below)  and  that 
outflow  boundary  conditions  of  various  kinds  be  employed,  including 
specification  of  v = 0 , specification  of  v_  = 0 , and  specification 
of  v ~ 0 • Most  of  the  results  reported  below  were  obtained  by 
using  v = 0 . 

The  y direction  was  resolved  by  using  eight  Fourier  modes  in 
a spectral  representation  (Gottlieb  and  Orszag,  1977).  This  repre- 
sentation ensures  periodicity  and  gives  an  adequate  representation  of 
the  initial  stages  of  nonlinear  growth  before  bursts  actually  appear. 

The  bursts  cannot  be  resolved  by  using  only  eight  Fourier  modes  because 
they  are  very  localized  in  y , so  our  calculations  must  be  considered 
inadequate  after  the  calculations  predict  bursts  to  occur.  Neverthe- 
less, our  predictions  of  the  location  of  transition  should  be  good. 
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We  resolve  the  boundary  layer  z direction  by  using  an  algebraic 
mapping  together  with  a 33-mode  Chebvshev  polynominal  expansion  (Gottlieb 
and  Orszag,  1977).  Rigid  boundary  conditions  are  applied  at  z = 0 , 
but  no  boundary  conditions  (other  than  boundedness)  are  applied  at 
z - 00  . The  advantages  of  Chebvshev  polynominal  solution  of  the 
z direction  are  that  it  is  highly  accurate,  relatively  easy  to  implement, 
and  gives  excellent  resolution  of  any  thin  boundary  layers  that  may 
appear. 

Time  stepping  is  done  by  a combination  of  methods.  A semi- implicit 
algorithm  is  used  in  x to  avoid  stringent  time-step  restrictions 
imposed  by  outflow  boundary  layers  (Gottlieb  and  Orszag,  1977,  Haidvogel 
and  Orszag,  1977).  The  semi-implicit  scheme  is  combined  with  Adams- 
Bashforth  explicit  second-order  time  differencing  on  the  advective 
terms  and  Crank-Nicolson  implicit  time  differencing  on  the  viscous 
terms.  The  overall  time-stepping  restrictions  are  not  too  severe  and 
are  of  the  form  At  ( Ax/u'  , where  Ax  is  the  spatial  resolution  in 
the  interior  of  the  computational  domain  and  u'  is  the  magnitude  of 
a typical  fluctuation  velocity.  In  practice,  accurate  numerical  results 
(as  opposed  to  just  stable  ones)  are  obtained  by  taking  the  time  step 
nearly  100  times  smaller  than  this  stability  bound.  To  resolve  accuratel 
the  advection  effects  of  the  unperturbed  mean  flow,  we  must  take  time 
steps  that  are  comparable  to  the  advection  time  of  the  mean  velocity  over 
one  effective  grid  interval. 

The  pressure  is  solved  by  fast  Fourier  transform  methods  in 
x and  v and  a matrix  inversion  method  in  z . The  matrix  multipli- 
cation that  is  required  is  fast  and  accurate;  the  pressure  computation 
requires  about  25%  of  the  overall  time  step. 

The  initial  conditions  we  have  used  are  those  described  bv  Orszag 
(1974).  They  consist  of  a superposition  of  a two-dimensional  Tollmien- 
Schlicting  wave  and  a three-dimensional  Tollmien-Schlicting  wave  super- 
posed on  a Blasius  profile,  as  suggested  by  the  laboratory  experiments  of 
Klebanoff,  Tidstrom  and  Sargent  (1962)  and  the  theory  of  Benney  and 
Lin  (Benney,  1964). 
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Transitlon  on  a Flat  Plate 

We  have  performed  a series  of  numerical  experiments  on  transition 

on  a flat  plate.  The  inlet  of  the  computational  domain  is  placed  1 m 

beyond  the  edge  of  a flat  plate  on  which  there  is  a zero-pressure  gradient 

free-stream  flow  incident  with  U = 1500  cm/s  . The  kinematic  viscosity 

2 

is  assumed  to  be  0.15  cm“/s,  the  boundary  layer  thickness  at  the  inlet 
is  6 = vx/U  = 0.099  cm,  the  Reynolds  number  at  the  inlet  is 
R = U6/v  = 1000  . The  length  of  the  computational  domain  is  chosen  to 
be  50  cm,  and  the  scale  height  in  the  algebraic  transformation  (2.1)  is 
chosen  to  be  H = 0.3  cm  . 

The  imposed  primary  Tollmien-Schlicting  wave  is  chosen  typically 
to  be  the  unstable  eigenmode  of  the  Orr-Sommerf eld  equation  with  wave- 
length 3.81  cm  and  has  an  imposed  maximum  u'  fluctuation  of  (0.01)U  . 
Various  three-dimensional  perturbations  on  this  primary  two-dimensional 
mode  are  investigated  below;  they  typically  consist  of  a (0.0015)U 
((15%)u')  perturbation  in  the  form  of  a three-dimensional  Tollmien- 
Schlicting  wave  with  the  same  x wavelength  as  the  primary  mode,  but 
with  various  spanwise  (y)  wavelengths. 

In  a typical  calculation,  the  time  step  is  chosen  to  be  (0.002)s, 
which  is  small  compared  with  the  primary  Tollmien-Schlicting  wave  period 
of  (0.007)s.  The  calculations  proceed  until  a statistically  steady 
state  develops,  normally  after  several  hundred  time  steps. 

In  figure  2 we  plot  the  results  of  a two-dimensional  flow  simulation 
with  only  the  two-dimensional  Tollmien-Schlicting  wave  applied.  Observe 
that  there  is  no  explosive  growth  of  the  amplitude  with  increasing 
x , which  indicates  absence  of  transition. 

In  figure  3 we  plot  the  results  of  a three-dimensional  simulation 
in  which  the  spanwise  wavelength  of  the  three-dimensional  perturbation 
is  2.5  cm.  In  figure  4 we  plot  the  results  of  a similar  experiment  in 
which  the  outflow  boundarv  conditions  u = v = w =0  are  replaced 

' XXX 

bv  u = 0 , and  v = w = 0 . 

x 

In  figures  5-7  we  compare  the  mean  velocity  profiles  at  three 
locations  in  the  flow  whose  perturbation  amplitudes  are  plotted  in 
figure  3.  Observe  that  as  the  downstream  distance  increases,  the  flow 
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becomes  highly  inflectional  at  the  spanwise  location  of  the  maximum 
velocity  fluctuation.  Beyond  the  point  at  which  these  highly  inflec- 
tional profiles  appear,  we  expect  the  rapid  occurrence  of  bursts. 

In  figure  8 we  show  the  effect  of  varying  the  spanwise  wavelength 
on  the  location  of  transition.  Observe  the  relatively  broad  minimum 
at  a spanwise  wavelength  of  2 cm.  Also,  observe  that  for  a spanwise 
wavelength  of  7 cm,  no  transition  was  observed  within  the  computational 
box. 

More  details  on  these  numerical  experiments  are  being  prepared  for 
publication  (Orszag,  1977).  Preliminary  results  have  been  published  in 
Orszag  (1976). 

Effect  of  Surface  Curvature 

We  have  investigated  the  effect  of  surface  curvature  on  the 
transition  process  by  solving  a modified  set  of  Navier-Stokes  equations 
that  include  the  primary  effects  of  such  a curvature.  Thus,  in  cylin- 
drical coordinates,  with  the  approximation  that  the  radius  r is 
constant  throughout  the  boundary  layer  (which  is  justified  if  r >>  5 , 
the  thickness  of  the  boundary  layer),  we  obtain 
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where  7“  also  includes  some  terms  that  depend  on  R , and  R is  the 
constant  radius  r . The  advantage  of  eqs.  (3. 1)— (3.4)  is  that  they  may 
be  solved  in  the  same  planar  geometry  used  for  the  flat  plat  transition 
experiments  reported  above.  The  boundary  conditions  are  still  u = v = w = 0 
at  the  flat  plate,  and  the  same  numerical  techniques  apply  to  the  solution 
of  the  problem.  The  only  difference  is  that  the  mean  inflow  generates 
longitudinal  vorticity  because  of  the  curvature  terms,  and  this  longitudi- 
nal vorticity  must  be  included  in  the  numerical  calculations.  Thus,  we 
ran  the  code  for  about  200  time  steps  to  allow  the  longitudinal  vorticity 
to  develop,  and  then  we  applied  the  Tollmien-Schlicting  waves  (two-  and 
three-dimensional)  as  perturbations. 


A 
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In  figure  9 we  show  the  results  of  a numerical  experiment  that  is 
otherwise  identical  to  the  experiment  plotted  in  figure  3,  except  that 
it  includes  the  effects  of  concave  surface  curvature  with  R = 15  cm  . 
Observe  that  the  location  of  transition  has  been  profoundly  affected  even 
though  R is  very  large  compared  with  the  boundary  layer  thickness  of 
0.1  cm.  Additional  results  of  this  kind  are  given  by  Orszag  (1977). 


uile  ol  the  maximum  perturbation  versus  x in  a calculation  o£  transition 
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4 . Galerkln  Approximation 

The  idea  of  the  Galerkin  approximation  method  discussed  here  is 
similar  to  that  of  nonlinear  stability  theory.  The  nonlinear  stability 
theory  of  boundary  layer  flows  was  worked  out  by  Benney  (1964)  and 

subjected  to  numerical  test  by  Antar  and  Collins  (1975).  Here  we  use  a 

slightly  different  approach.  ^ 

The  basic  idea  is  that  the  fundamental  nonlinear  dynamics  are  not 
grossly  affected  by  the  inflow-outflow  boundary  conditions  imposed  on 
the  full  simulations  discussed  in  section  3.  Therefore,  we  apply  the 

much  simpler  periodic  boundary  conditions  in  x and  y while  retaining 

the  rigid  boundary  conditions  imposed  in  the  z direction.  In  this  case, 
we  can  seek  the  solution  of  the  Navier-Stokes  equations  as  a Fourier 
series  in  x and  y as  follows: 

v (x , y , z , t)  = u(k,p,z,t)e  . (4.1) 

1 k | <K  1 p 1 <P 

Next,  the  Fourier  representation  (4.1)  is  substituted  into  the 
Navier-Stokes  equations,  and  equations  for  the  Fourier  components  are 
obtained  by  equating  coefficients  of  the  various  Fourier  modes.  Finally, 
a low-order  Galerkin  approximation  is  obtained  by  retaining  only  the 
lowest  order  modes  and  their  harmonics,  so  the  truncation  K = P = 2 is 
applied  to  (4.1). 

The  resulting  Galerkin  approximation  equations  are  solved  numerically 
as  differential  equations  in  time,  and  the  results  are  transformed  back 
to  a spatial  representation  by  using  a phase  velocity  transformation 
(even  though  the  group  velocity  transformation  is  more  justifiable). 

The  results  of  a comparision  with  the  three-dimensional  simulation  of 
figure  3 are  given  in  figure  10.  Although  the  details  of  the  explosive 
growth  of  perturbation  amplitude  are  not  given  accurately  by  the  Galerkin 
approximation,  the  qualitative  agreement  is  impressive.  In  particular, 
it  appears  that  the  Galerkin  approximation  can  predict  the  location  of 
transition  to  within  approximately  20%  even  though  the  Galerkin  approxi- 
mation requires  about  two  orders  of  magnitude  less  computer  time  than  the 
direct  solution  of  the  Navier-Stokes  equations. 

4 


son  of  the  predictions  of  a Galerkin  calculation  with  K = P = 2 and  the  full 
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The  apparent  success  of  the  Galerkin  approximation  employed  here 
is  a very  hopeful  development.  In  particular,  the  full  numerical  simu- 
lations of  transition  are  probably  too  complicated  and  costly  to  perform 
on  any  realistic  three-dimensional  body.  The  Galerkin  results,  however, 
suggest  that  the  boundary  conditions  are  not  too  important  in  detail,  so 
complicated  geometries  can  be  studied  by  isolating  small  regions  and 
piecing  the  results  together. 

A more  detailed  development  of  the  Galerkin  approximations  is  being 
prepared  for  publication. 


t 


I 
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